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ABSTRACT 

A previously developed local-global stiffness matrix methodology for the response of a composite 
half plane, arbitrarily layered with isotropic, orthotropic or monoclinic plies, to indentation by a rigid par- 
abolic punch is further extended to accommodate the presence of layers with complex eigenvalues (e.g., 
honeycomb or piezoelectric layers). First, a generalized plane deformation solution for the displacement 
field in an orthotropic layer or half plane characterized by complex eigenvalues is obtained using Fourier 
transforms. A local stiffness matrix in the transform domain is subsequently constructed for this class of 
layers and half planes, which is then assembled into a global stiffness matrix for the entire multilayered 
half plane by enforcing continuity conditions along the interfaces. Application of the mixed boundary 
condition on the top surface of the half plane indented by a rigid punch results in an integral equation for 
the unknown pressure in the contact region. The integral possesses a divergent kernel which is decom- 
posed into Cauchy-type and regular parts using the asymptotic properties of the local stiffness matrix and 
a relationship between Fourier and finite Hilbert transform of the contact pressure. The solution of the 
resulting singular integral equation is obtained using a collocation technique based on the properties of 
orthogonal polynomials developed by Erdogan and Gupta. Examples are presented that illustrate the 
important influence of low transverse properties of layers with complex eigenvalues, such as those exhi- 
bited by honeycomb, on the load versus contact length response and contact pressure distributions for half 
planes containing typical composite materials. 


INTRODUCTION 

In this paper we investigate the response of an arbitrarily laminated composite half plane, contain- 
ing layers whose displacement field is characterized by complex eigenvalues, indented by a well- 
lubricated rigid punch of a parabolic profile. A solution to this class of problems has been provided previ- 
ously for half planes consisting of isotropic, transversely isotropic, orthotropic or monoclinic layers 
characterized by real eigenvalues using a generalized plane deformation formulation, Fourier transforms 
and the local-global stiffness matrix approach (Pindera, 1991; Pindera and Lane, 1993; Binienda and Pin- 
dera, 1994). Displacement fields characterized by real eigenvalues occur in the majority of advanced 
continuously-reinforced composite materials in use today (Pagano, 1970). However, in the case of honey- 
comb layers oriented in the manner shown in Fig. 1, the solution of the Navier’s equations under the con- 
straint of generalized plane deformation in the x-z plane is characterized by complex eigenvalues. The 
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complex eigenvalues arise due to dramatic differences in the honeycomb’s elastic moduli as will be dis- 
cussed later. Herein, we present the solution to a class of composites characterized by complex eigen- 
values and subsequently incorporate it into the local-global stiffness matrix formulation of contact prob- 
lems involving arbitrarily layered half planes. Numerical examples are provided that illustrate the effect 
of low transverse properties of layers with complex eigenvalues, such as those exhibited by honeycomb, 
on the contact load versus contact length response and contact presure distributions for certain half plane 
configurations containing typical composite materials. Exact, analytical solutions to contact problems 
involving such sandwich configurations, which do not appear to be presently available, have 
technologically-significant applications in the aircraft, automotive and marine industry. We note that 
complex eigenvalue solutions also appear in the analysis of piezoelectric laminates, as shown by Heyliger 
(1994) and Heyliger and Brooks (1995) using the Fourier series representation of the displacement fields. 

As a first step, we construct a local stiffness matrix in the Fourier transform domain for layers with 
displacement fields characterized by complex eigenvalues by relating the displacement components at the 
layer’s top and bottom surfaces to the corresponding traction components. The local stiffness matrix for 
the corresponding half planes is recovered as a special case. Assembly of the local stiffness matrices of 
the individual layers comprising the arbitrarily laminated half plane into the global stiffness matrix in a 
particular fashion ensures satisfaction of displacement and traction continuity at the common interfaces 
and the external boundary conditions. Application of the mixed boundary condition on the slope of the 
surface displacement in the contact region and traction-free requirement elsewhere on the boundary pro- 
duces an integral equation for the determination of the contact pressure. As described elsewhere (cf., Pin- 
dera and Lane (1993)), the local/global stiffness matrix approach naturally facilitates decomposition of 
this integral equation into singular and regular parts that, in turn, can numerically be solved using the col- 
location technique outlined by Erdogan (1969) and Erdogan and Gupta (1972). The decomposition of the 
integral equation uses asymptotic properties of the local stiffness matrix and a relation between Fourier 
and Hilbert transforms of the contact pressure. 


ANALYSIS 

We consider a configuration composed of a number of layers bonded to a half plane wherein each 
region exhibits (transversely) isotropic, orthotropic, and/or monoclinic properties, Fig. 1 . A monoclinic 
ply is obtained by rotating a transversely isotropic, unidirectional ply through an angle 0 about its out-of- 
plane axis. A local x-y-z coordinate system is placed in the center of each layer such that the x and y 
axes lie in the lamination plane and the z axis is perpendicular to this plane. The layered medium is infin- 
ite in the x-y plane and the loading is such that the problem is plane in the x-z coordinate system. For 
the bottom half plane, the local coordinate system is placed at the bounding surface. The assemblage is 
indented by a rigid, frictionless punch of a parabolic profile and we are interested in the applied load as a 
function of the contact length and the resulting normal stress distribution in the contact region. 
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Due to the presence of monoclinic (i.e. off-axis) plies, generalized deformation formulation is 
employed with the displacement components taken in the form, 

w = w(jc,z), v = v(x,z), w = w(x,z) (1) 

Using the above displacement field, strain-displacement relations, and constitutive equations for a 
transversely isotropic or orthotropic ply with fibers rotated about the 2-axis, in the equilibrium equations, 
the Navier’s equations for a generic monoclinic ply become, 

Cn u ,xx + C55W tZZ + C 16V ju- + C 4 5 v tZZ + (C 13 + C 5 5)w xz = 0 

^16 w ,xx + ^45^,zz (^36 ^45) w ,xz = ^ (2) 

(C13 + C 55 )w JCZ + (C 36 + C 4 s)v X 2 + C5 5 W >xt + C33 W (22 = 0 

where the barred stiffness elements C l} are related to the unbarred stiffness elements C XJ in the principal 
material coordinate system by the familiar transformation equations. We note that coupling exists 
between all the displacement components for a monoclinic layer. The above equations can be specialized 
for an orthotropic or (transversely) isotropic layer by setting the stiffness elements Ci6, C36 and C 4 5 to 
zero, and by replacing the remaining Cy’s by C^’ s. In this case, the out-of-plane displacement component 
v(jc,z) becomes uncoupled from the in-plane displacement components u(x,z) and w(x,z). 

The solution of the equilibrium equations for each layer must satisfy the external surface mixed 
boundary conditions as well as the interfacial traction and displacement continuity conditions at each 
interface. The external surface mixed boundary conditions ensure that the normal traction component 
is zero outside the contact region Ijc I > c, while inside the contact region the vertical displacement 
w(x,z) conforms to the profile of the punch. The contact length 2 c is assumed to be given and the resul- 
tant load P is calculated by integrating the contact pressure distribution over 2 c, provided that remains 
compressive in the contact region. The formulation is thus valid provided that no separation takes place 
between the punch and the half plane during continuing loading (Urquhart and Pindera, 1994 ). These 
mixed boundary conditions are given in the form, 

w \,x( x >+h j/ 2 ) =f(x) for Ijc I <c and 0^ = 0 for \x \ >c ( 3 a) 

a xz ~ o yz = 0 for -00 < jc < +00 ( 3 b) 

while the interfacial traction and displacements continuity conditions are, 

u k (x,-h k / 2 ) = u k+i (x,h k+l / 2 ) o'l l (x,-h k /l) = o l iZ\x,h M /l) i=x,y,z ( 4 ) 

where u k = ( w k , u k , v k ), is the displacement vector for the kth ply. 
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The solution to eqns ( 2 ) subject to the boundary and continuity conditions given by eqns ( 3 ) and ( 4 ) 
is obtained using Fourier transforms defined, along with the inverse transform, by 


4 2k 

where «*(s,z) is the displacement vector for the kth ply. Application of the Fourier transform reduces the 
system of partial differential equations, eqns (2), to a system of ordinary differential equations in z with 
the transform variable s appearing as a parameter. For a monoclinic layer these equations are, 

C55"zz- S 2 C\\U +C45V (Z2 - J 2 C 16 V -«(Ci3+C 5 5)w z =0 

c 45« zz - s 2 C l6 u + C44V 22 - ^ 2 C 66 v - «(C 36 +C 45 )iv z = 0 

-is(C \i+C 55 )u z - w(C 3 6+C 4 5)v z + C 33 vv zz -5 2 C 3 5W =0 (6) 

As indicated previously, the corresponding equations for orthotropic or transversely isotropic layers are 
obtained by setting Cy = Cy, with C 16 = C 36 = C45 = 0. 

The solutions to the above equations are sought in the form vv(*s,z) = w 0 (s)e s ^ y u(s 9 z ) = wo(^)e s?tZ , 
and v(s>z) = v 0 (s)e sXz y where X is determined from a characteristic equation whose form depends on the 
layer’s material symmetry as discussed below. If the thicknesses of the layers comprising the half plane 
are finite, the exponential terms in the transform domain solutions are expressed in terms of hyperbolic 
functions to facilitate construction of the local stiffness matrix for a given layer. 

The eigenvalues X’s for monoclinic, orthotropic and isotropic layers or half planes are obtained 
from the following equations (upon substitution of the assumed displacement solutions into eqns (6)), 

monoclinic : -A X 6 + BX 4 + CX 2 + D = 0 

orthotropic : (C44X 2 - C^XAX 4 + BX 2 + C) = 0 

isotropic : (X 2 - 1 )(X 4 - 2X 2 + 1) = 0 (7) 

where the coefficients A through D for are lengthy algebraic expressions involving the elastic stiffness 
elements Cy’s for monoclinic layers and Cy’s for orthotropic layers (cf. Pagano ( 1970 )). The eigenvalues 
for an isotropic layer are real (i.e., Xi 2 = ± 1 , X 3 4 = ±1, X 5 6 = ± 1 ). In the case of most advanced unidirec- 
tional composites modeled as either transversely isotropic or orthotropic, the eigenvalues X’s are also real. 
The same holds true for an orthotropic unidirectional composite laminae rotated through an off-axis angle 
in its plane (i.e., the x-y plane in Fig. 1 ). Such a rotated laminae behaves like a monoclinic layer in the 
fixed coordinate system x—y—z shown in Fig. 1 . Solutions for isotropic, orthotropic and monoclinic lam- 
inae with real eigenvalues, and subsequent construction of the local stiffness matrix, have been provided 
by Pindera and Lane ( 1993 ). In the case of orthotropic materials with honeycomb-type microstructures, 


j u k (s,z)e ,sx ds ( 5 ) 


1 


+00 


u k (s y z) = -prr f U k (x>z)e tsx dx , u k (x y z) 
v2tc jL 
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however, the eigenvalues are typically complex. The solution for this class of materials is provided next. 


Solution for orthotropic materials with complex eigenvalues 

The coefficients A , fi, and C for an orthotropic material appearing in the second of eqn (7) are given 
in terms of the elastic stiffness constants as follows, 

A = C33C55 , B = C 13 (C 13 + 2C 55 ) - C\\ C33 , C = C n C 55 


The expressions for the eigenvalues of an orthotropic material are thus obtained in the following form: 



C\\ C33 - C 13 (C 13 + 2C55) + yj[C ] . C33 - C 13 (C 13 + 2C55)] 2 *“ 4 C n C33C55 

2C33C55 



^11^33 ~ c 13(^13 + 2 c 55 ) - V[c 1 1 C33 - c i3(c i3 +2C55)] 2 -4C11C33C55 

2C33C55 


^5,6 “ ±VC 6 6/C44 (8) 

Solutions with complex eigenvalues to eqn (6) are obtained when the expression under the square root 
appearing in the first two sets of eigenvalues in eqn (8) becomes negative. This expression can be fac- 
tored to determine the region in the elastic stiffness constant space where complex eigenvalues are found. 
Setting this expression to zero, 

(C 13 + VC 1 1 C33 )(C 13 - VC11C33 )[C [3 - (-2 C55 4* VC11C33 )][C 13 - (-2C55 - V^i 1 C 33 )] = 0 

yields four planes in the ^jC {] C 33 -C 13 -C 5 5 space that separate the region with real and complex eigen- 
values, Fig. 2. These planes are: C 13 = ±^jC\[C^ , C 13 = -2C55 ± VC11C33 . 

Complex expressions for the eigenvalues in the solutions for the displacement components u(s,z) 
and w(s,z) are obtained for honeycomb materials which are characterized by very large Young’s modulus 
associated with the z direction, £33, and small Young’s moduli associated with the x and y directions, £n 
and £22 1 respectively. In this case the eigenvalues have the form, 

2 = a ±ib , ^3 4 ~~a ± ib (9) 

where 

a = Vrcos(j) , b = Vr sin(y) , r = 

0 . V4 CuCnCls - [C,,C 3 3 -C I3 (C 13 +2C 55 )] 2 

0 = tan ( 

C11C33 - C 13 (C,3 +2C55) 


) 


( 10 ) 



Thus the complex eigenvalue solutions to the Fourier-transformed Navier’s equations for orthotro- 
pic layers are conveniently expressed in terms of hyberbolic functions as follows, 

u(s,z) = [ Fi(s)cos(sbz ) + G i(s)sin (sbz)]sinh (saz) + [F 2 (s)cos (sbz) + G 2 (s)sin(sbz)]cosh(saz ) 

= i [/ ! (s)cos (sbz) + K { (s)sin (sbz)]sinh (saz) + [/ 2 (.r)coj (sbz) + K 2 (s)sin (sbz)]cosh (saz) 

v(s,z) = F 2 (s)cosh(sX^z) + G^(s)sinh(sXsz) (1 1) 

where the coefficients l,(s) and K,(s) are related to the coefficients F t (s) and Gj(s) (i = 1, 2) as follows, 

I i( s )=fo^2( s ) + So^i ( s ) » h( s ) = foFi( s ) + 8oG 2 (s) 

K\(s) — -go^iis) +foG 2 (s ) , K [(s) = -g 0 F 2 (s) + f 0 G i(s) (12) 

where 


a[(a 2 +fc 2 )C55-C n ] fc[(q 2 +fc 2 )C 5 5 + C, 1 ] 

(a 2 + b z )(C + C 55 ) (q 2 +fo 2 )(C| 3 + C 55 ) 

For orthotropic half planes, on the other hand, it is convenient to express the corresponding solu- 
tions in terms of exponential functions as follows, 

u(s,z)= [/ r |(s)c0j(l.y I bz) +F 2 (s)sin(\s I bz)]e u ' az 

w(s,z) = sgn(s)i[(f 0 F i(s) + g 0 F 2 (s))cos ( Is \bz) + (-goF^s) +f 0 F 2 (s))sin ( I s I bz)]e luz 

v(s,z) = F 3 Cy)e IXjZ (14) 

The displacements given by eqn (11) are substituted back into the constitutive equations in order to 
determine interfacial tractions needed in applying the continuity conditions, eqn (4), in the transform 
domain, and ultimately in the construction of the local stiffness matrix. For an orthotropic layer, the con- 
tinuous interfacial stresses given in terms of displacements in the transform domain are, 

O a (j,z) = -wC|3M(^z) + C 3 3VV jZ (s,z) 

O xz (s,z) = c 55 (u (s,z) - isw(s,z)) 

Gyz(s>Z) = C U V, z ( 15 ) 

Rather than determining the unknown Fourier coefficients Fj(s) and Gj(s) in the solutions for the dis- 
placement field in terms of the unknown contact pressure distribution through the application of the 
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interfacial continuity conditions and the mixed boundary condition at the top surface, the problem is 
reformulated in terms of the unknown interfacial displacements using the local-global stiffness matrix 
approach outlined in the sequel. 

Local-Global Stiffness Matrix Formulation 

The local-global global stiffness matrix approach naturally facilitates the reduction of the contact 
problem to a singular integral equation in the unknown pressure distribution in the contact region. This 
singular integral equation is then solved using the technique developed by Erdogan (1969) and Erdogan 
and Gupta (1972). First, a local stiffness matrix is developed for the Icih layer that relates the traction 
components on top and bottom surfaces of the layer to the corresponding displacement components. This 
is accomplished by expressing the Fourier coefficients Fj(s) and Gj(s) in terms of interfacial displace- 
ments obtained from eqns (11), and subsequently using these expressions in the interfacial tractions 
obtained from eqns (14). For orthotropic layers with complex eigenvalues, the resulting local stiffness 
matrix equation has the following form, 
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where the 'V and superscripts refer to the top and bottom surfaces of a given ply, respectively. The 
above equation is expressed symbolically in the form. 



^12 


f " 

ut 
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K 22 

4 

v~k 

► zz’i 

T k 

L 



< * 


„ j 


( 17 ) 


The elements of the local stiffness matrix are functions of the transform variable 5 , material stiff- 
ness constants C /y , and layer’s geometry. These are provided in the Appendix. The elements k^ for 
monoclinic, orthotropic, and transversely isotropic layers characterized by real eigenvalues have been 
provided previously by Pindera (1991) and Pindera and Lane (1993). 

Next, imposition of continuity of displacements and tractions along the common interfaces together 
with the external boundary conditions (eqns (3)-(4)) gives rise to a system of equations in the unknown 
interfacial displacement components represented in the matrix form below. It is observed that the assem- 
bly of the global stiffness matrix for the entire layered medium is carried out by superposing local 
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stiffness matrices of the individual layers along the main diagonal of the global matrix in an overlapping 
fashion. 


*11 

*12 

*2. 

K\ 2 + K 

0 

* 2 . 

0 

0 

0 




*22 +*11 


k) 


21 



( 18 ) 


In the above, K\{ is the stiffness matrix of the bottom half plane. 

Inverting the global stiffness matrix yields a relation between the top surface displacements and the 
top surface tractions, 



= [HUU T x 


( 19 ) 


where H\\ is the first 3x3 submatrix of the inverse of the global stiffness matrix in eqn (18). Defining 
the elements of the matrix H \ j by Hy, the transform of the normal displacement on the surface of the half 
plane can be expressed in terms of the normal contact stress in the absence of friction (i.e. 

<5xz — Gy Z — 0), 

1 1 1 +c 

w j =H {l (s)o zz (s)/s with G zz {s)=p(s) = -=[ p(x')e lsx dx 

"V —c 

Imposing the top surface mixed boundary condition on the slope of the normal displacement, w \ x - f(x ) 
in the interval -c < x < c, an integral equation for the unknown contact stress distribution p(x) is then 
obtained in the form, 


oo > oo 

w u = -=L J JW, e~ hx ds = j-j H u (s)p(s) e~ isx ds (20) 

The above is a singular integral equation because as s approaches ± 00 the kernel H\\(s) does not vanish, 
making the integral unbounded. To find the unknown function p(x), the dominant or singular part of the 
kernel must be identified. This is accomplished by first examining the asymptotic behavior of the local 
stiffness matrix as the transform variable s approaches positive or negative infinity. In this case, the cou- 
pling submatrices K\ 2 and in eqn (17) vanish and the stiffness matrix assumes the following form, 
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( 21 ) 
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K \\ 
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ut 
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The elements of the asymptotic stiffness matrices K\\ and K 22 are functions of the material properties of 
the given layer but not functions of the transform variable s or geometry (Pindera (1991)). 

The asymptotic behavior of the local stiffness matrix given by eqn (21) ensures that in the limit as s 
approaches ± °° , the resulting global stiffness matrix has only diagonal elements, so that there is no cou- 
pling between top and bottom surfaces of each lamina of the layered half plane for this limiting case. 
Consequently, the limiting behavior of eqn (19) becomes 

|(/ij=[^i 1 r 1 |fij (22) 


These results are used to separate the divergent integral in eqn (20) into one integral containing a 
Cauchy kernel and another integral with a regular kernel. Defining sgn(s) H* u = lim H u (s) to be the 

X— >±P° 

first element of the inverse of the asymptotic stiffness matrix in eqn (22) (or eqn (19)), where //,,(*) is 
the first element of the inverse of the global stiffness matrix in eqn (18), the mixed boundary condition 
given by eqn (20) becomes, 


Wl,x= 2n f sgn(s)Hu p(s) e ux ds - j ( H ,,(*) - sgn(s) H* u )p(s) e~ isx ds 


2k 


(23) 


In view of eqn (22), a relationship between the Fourier and finite Hilbert transforms of the contact pres- 
sure can be derived in the following form (cf. Gladwell (1980, p.210), 

—= J - i sgn (s) p(s ) e~ lsx ds = — f dt (24) 

<2l t_oo Ki c t-x 

reducing the dominant term of the singular integral to an integral containing a Cauchy kernel. Using odd- 
even properties of the regular kernel, the following form of the singular integral equation for the contact 
stress distribution is then obtained, 

H* +C /v 00 +C 

w l.x = '77 = — 7~j “■<* + “/ J H°n(s)p(t) sin(r-jc)j dt ds (25) 

K K _ £ t X 71 o -c 

where R is the punch radius and //?| (s) - H j [ ( 5 ) - sgn{s)H\\ is the regular kernel. 
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NUMERICAL RESULTS 

The solution to eqn (25) has been obtained using the numerical technique for singular integral equa- 
tions developed by Erdogan and Gupta (1969, 1972) which is based on orthogonal properties of Che- 
byshev polynomials in a Gaussian integration approach. The details of the application of the technique to 
the outlined contact problem have been provided by Pindera and Lane (1993). 

Here, we illustrate the developed solution by investigating the effect of stacking sequence and 
material properties on the contact load as a function of contact length, and normalized contact stress pro- 
file, for layered half planes constructed with commonly used metal matrix and polymeric matrix unidirec- 
tional composites, and honeycomb-like layers. The unidirectional composites employed in the analysis 
include three types of polymeric matrix composites, namely glass/epoxy (Gl/Ep) and two types of 
graphite/epoxy (Gr/Ep) (T300/934 or Gr/Ep-1, and P75/934 or Gr/Ep-2), and two types of metal matrix 
composites, namely boron/aluminum (B/Al) and P100/6061 graphite/aluminum (Gr/Al). These compo- 
sites provide a wide range of material properties with different orthotropy ratios. The elastic constants of 
these materials are given in Table L The honeycomb properties included in the table are based on the data 
reported by Shuart (1978) for an aluminum honeycomb whose macroscopic elastic properties were gen- 
erated using a finite-element homogenization analysis. For the purpose of the present illustration, these 
properties (except the Poisson’s ratios) were multiplied by a factor of seven in order to avoid excessively 
high deformations in the contact region caused by the low value of the Young’s modulus in the direction 
of the applied load, which may have invalidated the present analysis based on linear elasticity. We point 
out that the Young’s modulus in the direction of the applied load of the stiffer honeycomb is now of the 
same order of magnitude as that of glass/epoxy, graphite/epoxy and graphite/aluminum, and an order of 
magnitude smaller than that of boron/aluminum. We also note that, in practice, honeycomb properties can 
be controlled by the choice of material, wall thickness and cell dimensions. Figure 3 shows the locations 
of three of the six composite materials and the honeycomb in the VC 11 C 33 -Ci 3 -C 55 space, with the 
planes that define the region with complex eigenvalues included. 

For the layered half planes, the following configurations were studied: a honeycomb layer bonded 
to a composite half plane; a composite layer bonded to a honeycomb half plane; and a composite layer 
bonded to a honeycomb layer which, in turn, was bonded to a composite half plane. We note that 
although the configuration with the honeycomb layer on top is neither practical nor technologically 
important, nevertheless the results can be used for comparison with more technologically meaningful 
configurations as well as correlation with past results. The composite layers in all the investigated confi- 
gurations had fibers oriented along the Jt-axis (see Fig. 1). Thus their response in the x-z plane was that of 
an orthotropic material characterized by real roots in the displacement field solutions. The thickness of 
the surface and interior plies comprising the layered half planes was 1 .27 and 2.54 mm, respectively. The 
punch radius employed in the calculations was 25.4 mm. 

As a first step, and to provide a point of reference, we compare the load versus contact length 
responses of homogeneous half planes constructed with the aforementioned materials, Fig. 4. As dis- 
cussed previously by Binienda and Pindera (1994), the response of the homogeneous composite half 
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planes with fibers oriented along the x-axis is primarily influenced by the Young’s modulus in the direc- 
tion of the applied load, £ 33, with the longitudinal modulus E\\ playing a less significant role. Thus the 
contact load versus contact length response is stiffest for the composite with the largest £33, and 
decreases with decreasing £33. This trend, however, does not hold for the honeycomb half plane which 
exhibits the most compliant response despite its larger £33 relative to that of the two Gr/Ep half planes. 
Evidently, the properties transverse to the load direction play a substantial role in this case. To 
emphasize this point, we have included in Fig. 4 the response of a fictitious isotropic material, which will 
be referred to as "isotropic honeycomb", with the Young’s modulus equal to the honeycomb’s modulus in 
the direction of the applied load, £33, and a Poisson’s ratio of 0.30. Considerable stiffening of the load 
versus contact length response is observed by increasing the Young’s moduli transverse to the load direc- 
tion to equal that in the direction of the applied load. 

Next, we compare the response of the configuration comprised of the honeycomb layer bonded to a 
composite half plane, Fig. 5, to that of the configuration comprised of a composite layer bonded to the 
honeycomb half plane, Fig. 6. In the first instance, the initial load versus contact length response, Fig. 5a, 
follows closely the response of the honeycomb half plane irrespective of the properties of the supporting 
half plane. Thus, as is well known, the initial response is governed by the properties of the top layer. 
Beyond the contact half length of approximately 1 mm, the influence of the supporting half plane is 
becoming apparent, with the trends following those seen in Fig. 4. The normalized pressure distributions 
for these configurations at the contact half length of 2.54 mm, i.e., c/h j = 2 (where h \ is the thickness of 
the top layer, see Fig. 1), given in Fig. 5b, exhibit departures from elliptical profiles that increase with 
increasing Young’s modulus in the direction of the applied load. These departures are characterized by a 
higher maximum pressure in the center of the contact region (i.e., x/c = 0.0) relative to that of an elliptical 
profile. When the stacking sequence is reversed in the second instance. Fig. 6, the influence of the com- 
posite surface layer’s properties is immediately evident in the load versus contact length response, Fig. 
6a. The configuration with the stiffest surface layer and lowest £n/£ 3 3 ratio, i.e., B/AI-honeycomb half 
plane, exhibits load versus contact length response that initially departs substantially from parabolic, indi- 
cating the presence of localized bending of the surface layer. This is also evident in the contact pressure 
distribution, Fig. 6b, which is pronouncedly nonelliptical, characterized by maximum values occuring at 
the edges of the contact region. The remaining configurations exhibit parabolic load versus contact length 
responses and nearly elliptical contact pressure distributions. As in the previous case, the small departures 
from elliptical profiles increase with increasing Young’s modulus in the direction of the load and decreas- 
ing ratios £n/£33. In this case however, the departures are characterized by a lower maximum pressure 
at the center of the contact region relative to that of an elliptical profile. This, in turn, indicates decreasing 
resistance to localized bending, with the B/Al-honeycomb half plane being the least resistant. 

The substantial impact of the low transverse properties of the honeycomb substrate on the contact 
responses shown in Fig. 6 is highlighted by comparison with the corresponding results, given in Fig. 7, 
generated using the properties of the fictitious isotropic honeycomb introduced earlier. Comparison of 
Figs. 6a and 7a demonstrates the considerable stiffening of the contact load versus contact length 
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response due to the increase in the transverse properties of the honeycomb substrate. The extent of the 
stiffening is most evident for the B/Al configuration which nearly regains a parabolic contact load versus 
contact length response due to the increase in the transverse properties of the substrate. Comparison of 
Figs. 6b and 7b demonstrates that the departures from elliptical contact pressure distributions are reduced 
at the given contact length (c = 2.54 mm) due to this increase, which is particularly evident for the B/AL 
configuration. 

The last investigated configuration is composed of a composite layer bonded to the honeycomb 
layer which, in turn, is bonded to a composite half plane. Figures 8a and 8b present the load versus con- 
tact length response and the contact pressure distributions at the contact half length of 2.54 mm, respec- 
tively, for these configurations. The load versus contact length responses appear parabolic, and are stiffer 
than the corresponding responses observed in the configurations composed of a composite layer bonded 
to the honeycomb half plane. The influence of the supporting half plane is thus very much in evidence in 
these configurations. The contact pressure distributions are nearly elliptical for most configurations with 
the exception of the configuration with the B/Al ply. In this case, the pressure distribution is similar to, 
but not as pronouncedly nonelliptical as, that seen in Fig. 6b, suggesting the presence of some localized 
bending despite the parabolic load versus contact length response. Replacing the honeycomb layer with 
the fictitious isotropic honeycomb layer (not shown) does not visibly stiffen the contact load versus con- 
tact length response and decreases the departures from elliptical contact pressure distributions at the con- 
sidered contact length only for the B/Al configuration, in contrast with the preceding configuration. This 
demonstrates that the extent of the influence of honeycomb layers’ low transverse properties on the con- 
tact response of layered half planes also depends on the layer dimensions. 


CONCLUSIONS 

The capability to analyze arbitrarily layered half planes with differently oriented composite plies, 
indented by a frictionless, rigid parabolic punch, was extended to enable incorporation of layers or half 
planes characterized by complex-eigenvalue displacement fields into the analysis. Honeycomb layers or 
half planes with very low elastic properties perpendicular to the applied load fall into this category of 
materials. Expresssions for the elements of the local stiffness for such materials were developed and 
incorporated into the solution strategy for the contact problem of arbitrarily layered half planes based on 
the local/global stiffness matrix approach in the Fourier-transform domain. 

Unlike homogeneous half planes constructed with typical advanced unidirectional composites, the 
load versus contact length response of honeycomb-like homogeneous half planes is significantly influ- 
enced by the elastic constants associated with directions perpendicular to the applied load. The low values 
of these constants relative to the Young’s modulus in the direction of the applied load, E 33 , i.e., along the 
axis of the honeycomb, substantially degrade the load versus contact length response relative to that of 
homogeneous composite half planes with comparable Young’s moduli £ 33 but higher £,,. When a 
honeycomb-type layer is inserted directly underneath the top layer of a half plane laminated with typical 
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advanced unidirectional composites, the ability of the top layer to resist localized bending under the 
punch tends to be degraded due to the honeycomb’s low transverse properties. This may result in nonpar- 
abolic contact load versus contact length responses and nonelliptic pressure distributions with maximum 
magnitudes occuring at the outer edges for surface layers with sufficiently high values of £ 33 and 
E \\/E 33 . Increasing the transverse Young’s moduli of the honeycomb-like layer to equal the modulus in 
the direction of the applied load was shown to substantially stiffen the load versus contact length 
response, and decrease the departure of the pressure distributions from elliptical, thereby directly demon- 
strating the importance of honeycomb layers’ transverse properties. 

The presented solution methodology can be employed in studying the contact response of laminated 
plates containing layers with complex eigenvalues bonded to a stiff foundation whose response can be 
approximated by a half plane. Symmetrically laminated finite-thickness plates loaded by opposing sur- 
face contact loads, such as those exerted by aligned rollers, can also be investigated by modifying the 
boundary conditions through the imposition of zero vertical displacement along the plane of plate’s sym- 
metry. Similary, finite-thickness laminated plates supported at the bottom surface by a system of forces 
can be investigated by modifying the formulation to include the effect of the supports in the manner 
described by Pindera and Lane (1993). In this case, the singular integral equation for the contact stress 
distribution contains the unknown support reactions, requiring an iterative solution approach. 

The contact pressure distribution is the first step in obtaining a solution to layered half planes and 
plates indented by a punch. The knowledge of this distribution allows the calculation of internal stresses 
to determine the regions where potential damage may occur due to the imposition of concentrated loads. 
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APPENDIX - Elements of Local Stiffness Matrix 


The elements of the local stiffness matrix for orthotropic layers whose generalized plane deformation solutions in 
the Fourier transform domain are characterized by complex eigenvalues are given below. 


1 L I ^2 

ku=--CnP x {— + — ) 

2 Aj A 2 

, 1 r (“Cl 3 + C 33 P 3 L 3 ) (-C| 3 +C 33 P 3 L 4 ) 
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The asymptotic expressions for the elements of the local stiffness matrix for orthotropic layers with complex eigen- 
values as the transform variable s — » and for the corresponding half planes are given below. 


±&n -±C^(a-\ b ) 
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33 
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So 
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where the notation ± K*j denotes limiting behavior of kjj as 5 goes to ± °° 
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Table 1. Material properties of polymeric and metal matrix unidirectional composites. 


Material property 

Gl/Ep 

Gr/Ep- 1 
(T300/934) 

Gr/Ep-2 

(P75/934) 

B/Al 

Gr/Al 

Honeycomb 

E x x (GPa) 

42.7 

144.8 

243.0 

227.5 

402.6 

1.2 xlO ' 2 

E 22 (GPa) 

11.7 

10.3 

7.2 

137.9 

24.1 

0.7 xlO ' 2 

£33 (GPa) 

11.7 

10.3 

7.2 

137.9 

24.1 

11.4 

V 12 

0.27 

0.30 

0.33 

0.24 

0.29 

1.10 

V 13 

0.27 

0.30 

0.33 

0.24 

0.29 

0.35 xlO ' 3 

v 23 

0.55 

0.50 

0.49 

0.40 

0.45 

0.21 xlO ' 3 

G 12 (GPa) 

8.24 

5.51 

3.93 

55.15 

16.75 

1.52 

G ,3 (GPa) 

8.24 

5.51 

3.93 

55.15 

16.75 

4.66 

G23 (GPa) 

3.78 

3.45 

2.41 

49.24 

8.34 

1.95 
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Figure 1. Multilayered composite half plane with honeycomb layers by a 
rigid, parabolic punch. 
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Figure 2. Planes separating regions with real and complex eigenvalues of 
eqn. (8) in the ^/C n C 33 - C 13 - C 55 space. 
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Figure 3. Locations of the values of the stiffness elements Cy of several 
common unidirectional composites and an aluminum honeycomb in the 
•n/C m C 3 3 - C l3 - C 55 space. 
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Figure 5. Contact response of an aluminum honeycomb layer bonded to a 
composite half plane: a) load vs contact length; b) contact pressure 
distribution for c=2.54 mm. 
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Figure 7. Contact response of a composite layer bonded to an isotropic 
half plane with the same Young’s modulus as the honeycomb’s modulus in 
the direction of the applied load: a) load vs contact length; b) contact 
pressure distribution for c=2.54 mm. 
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Figure 8. Contact response of composite and aluminum honeycomb layers 
bonded to a composite half plane: a) load vs contact length; b) contact 
pressure distribution for c=2.54 mm. 


25 





REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of Information. Send comments regarding this burden estimate or any other aspect of this 
collection of Information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. 


1. AGENCY USE ONLY (Leave blank) 


4. TITLE AND SUBTITLE 


REPORT DATE 

December 1997 


3. REPORT TYPE AND DATES COVERED 


Frictionless Contact of Multilayered Composite Half Planes Containing Layers 
With Complex Eigenvalues 


AUTHOR(S) 

Wang Zhang, Wieslaw K. Binienda, and Marek-Jerzy Pindera 


Final Contractor Report 


5. FUNDING NUMBERS 


WU-523-2 1-1 3-00 
NAG3-1827 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESSEES) 

University of Akron 

Civil Engineering Department 

Akron, Ohio 44325 


0. PERFORMING ORGANIZATION 
REPORT NUMBER 


E-11008 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


NASA CR— 97-206309 


11. SUPPLEMENTARY NOTES 


Project Manager, Steven M. Arnold, Structures and Acoustics Division, NASA Lewis Research Center, organization 
code 5920, (216) 433-3334. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category: 24 


12b. DISTRIBUTION CODE 


Distribution: Nonstandard 


This publication is available from the NASA Center for AeroSpace Information, (301) 621-0390. 


13. ABSTRACT (Maximum 200 words) 

A previously developed local-global stiffness matrix methodology for the response of a composite half plane, arbitrarily 
layered with isotropic, orthotropic or monoclinic plies, to indentation by a rigid parabolic punch is further extended to 
accommodate the presence of layers with complex eigenvalues (e.g., honeycomb or piezoelectric layers). First, a general- 
ized plane deformation solution for the displacement field in an orthotropic layer or half plane characterized by complex 
eigenvalues is obtained using Fourier transforms. A local stiffness matrix in the transform domain is subsequently con- 
structed for this class of layers and half planes, which is then assembled into a global stiffness matrix for the entire 
multilayered half plane by enforcing continuity conditions along the interfaces. Application of the mixed boundary 
condition on the top surface of the half plane indented by a rigid punch results in an integral equation for the unknown 
pressure in the contact region. The integral possesses a divergent kernel which is decomposed into Cauchy-type and 
regular parts using the asymptotic properties of the local stiffness matrix and a relationship between Fourier and finite 
Hilbert transform of the contact pressure. The solution of the resulting singular integral equation is obtained using a 
collocation technique based on the properties of orthogonal polynomials developed by Erdogan and Gupta. Examples are 
presented that illustrate the important influence of low transverse properties of layers with complex eigenvalues, such as 
those exhibited by honeycomb, on the load versus contact length response and contact pressure distributions for half 
planes containing typical composite materials. 


14. SUBJECT TERMS 15. NUMBER OF PAGES 

Honeycomb composite; Contact problem; Frictionless; Indentation; Local/global stiffness; 31 

Singular integrals |i6. PRICE CODE 


17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 
OF REPORT OF THIS PAGE 


Unclassified 


Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


20. LIMITATION OF ABSTRACT 


Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. Z39-18 
298-102 






















